Assignment7_VIS

library(tidyverse)
library(sf)
library(leaflet)
library(htmltools)
library(htmlwidgets)
library(raster)
library(gstat)
library(spatial)
library(dplyr)
library(wesanderson)

leaflet_plane <- "+proj=longlat +datum=WGS84"

MA_state_plane <- "+proj=lcc +lat_1=41.71666666666667 +lat_2=42.68333333333333 +lat_0=41 +lon_0=-71.5 +x_0=200000 +y_0=750000 +ellps=GRS80 +units=m +no_defs"

traffic <- st_read("traffic.csv", quiet = T) %>%
  filter(Lat != "NA")

traffic2 <- st_as_sf(traffic, coords = c("Long", "Lat"), crs = crs(leaflet_plane)) %>%
  mutate(total_peak = as.numeric(as.character(Peak_total)))

streets <- st_read("http://bostonopendata-boston.opendata.arcgis.com/datasets/cfd1740c2e4b49389f47a9ce2dd236cc_8.geojson?outSR=%7B%22latestWkid%22%3A2249%2C%22wkid%22%3A102686%7D", 
                   quiet = T) %>%
  st_transform(crs = leaflet_plane)

south_boston <- st_read("C:/Users/mwill/OneDrive - Harvard University/2020 - Fall/Spatial Analysis/InteractiveMaps/SouthBoston_tract", quiet = T) %>%
  st_transform(crs = leaflet_plane)

leaflet(south_boston) %>%
  addProviderTiles(providers$Stamen.Toner) %>%
  addPolygons() %>%
  addPolylines(data = streets)

south_boston_streets <- streets[south_boston,]

## although coordinates are longitude/latitude, st_intersects assumes that they are planar

ggplot(south_boston_streets) +
  geom_sf()

south_boston_traffic <- st_join(south_boston, traffic2, join = st_intersects)

south_boston_traffic2 <- south_boston_traffic %>%
  group_by(GEOID, NAMELSAD, Associated.Streets) %>%
  summarise(total_peak_tract = sum(total_peak))

traf_mean <- mean(south_boston_traffic2$total_peak_tract, na.rm = T)

south_boston_traffic2$total <- ifelse(is.na(south_boston_traffic2$total_peak_tract) == T,
                                      traf_mean,
                                      south_boston_traffic2$total_peak_tract)

south_boston_traffic2$describe <- paste("<p>",south_boston_traffic2$NAMELSAD,
                                        "</p>Intersection:", 
                                        str_to_title(
                                          south_boston_traffic2$Associated.Streets),
                                        "<br>",
                                        "Vehicles during Peak Hours:",
                                        prettyNum(south_boston_traffic2$total_peak_tract,
                                                  big.mark = ",")) %>%
  lapply(htmltools::HTML)

bins <- seq(min(south_boston_traffic2$total_peak_tract, na.rm = T),
            max(south_boston_traffic2$total_peak_tract, na.rm = T), by = 1)
pal <- colorNumeric(palette = c("#fee6ce", "#fdae6b", "#e6550d"), 
                    domain = south_boston_traffic2$total_peak_tract)

leaflet(south_boston_traffic2) %>%
  addProviderTiles(providers$Stamen.TonerLite) %>%
  addPolygons(fillColor = pal(south_boston_traffic2$total_peak_tract), 
              fillOpacity = 0.4, 
              stroke = F, 
              highlightOptions = highlightOptions(fillOpacity = 1),
              popup = ~describe) %>%
  addLegend(pal = pal,
            values = ~total,
            title = "Total Vehicles",
            bins = 5,
            position = "topright")

pal <- colorNumeric(palette = c("#f2f0f7",
                                "#cbc9e2",
                                "#9e9ac8",
                                "#756bb1",
                                "#54278f"),
                    domain = traffic2$total_peak)
bins <- seq(min(traffic2$total_peak, na.rm = T),
            max(traffic2$total_peak, na.rm = T), by = 1000)

traffic2$describe <- paste("Intersection:", 
                                        str_to_title(traffic2$Associated.Streets),
                                        "<br>",
                                        "Vehicles during Peak Hours:",
                                        prettyNum(traffic2$total_peak,
                                                  big.mark = ",")) %>%
  lapply(htmltools::HTML)

leaflet(traffic2) %>%
  addProviderTiles(providers$Stamen.Toner) %>%
  addCircles(fillColor = pal(traffic2$total_peak), fillOpacity = .8,
             stroke = F, radius = 50,
             popup = ~describe,
             highlightOptions = highlightOptions(fillOpacity = 1)) %>%
  addLegend(pal = pal,
            values = ~total_peak,
            bins = 5,
            position = "topright",
            title = "Vehicles Counted during<br>Peak Hours 11/21/19")  

traffic_pts_sp <- traffic2 %>%
  st_transform(MA_state_plane) %>%
  as_Spatial()

tract_poly_sp <- south_boston %>%
  st_transform(MA_state_plane) %>%
  as_Spatial()

streets_poly_sp <- south_boston_streets %>%
  st_transform(MA_state_plane) %>%
  as_Spatial()

southie_raster <- raster(tract_poly_sp, res = 10)

gs <- gstat(formula = Peak_total~1, locations=traffic_pts_sp)
idw_interp <- interpolate(southie_raster, gs)

## [inverse distance weighted interpolation]

idw_interp_clip <- mask(idw_interp, streets_poly_sp)

bins <- seq(min(traffic2$total_peak),
            max(traffic2$total_peak), by = 1)
pal <- colorNumeric("viridis", 
                    domain = traffic2$total_peak,
                    na.color = "#00000000")

leaflet(traffic2) %>%
  addProviderTiles(providers$Stamen.Toner) %>%
  addRasterImage(idw_interp_clip, colors = pal, opacity = 0.8) %>%
  addLegend(pal = pal,
            values = ~total_peak,
            bins = 5,
            position = "topright",
            title = "Total Vehicles Counted<br>Peak Hours on 11/21/19")